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Abstract A viscoelastic model of the k-bkz (Kaye 1962; Bernstein et al. 1963) 
type is developed for isotropic biological tissues, and applied to the fat pad of the 
human heel. To facilitate this pursuit, a class of elastic solids is introduced through 
a novel strain-energy function whose elements possess strong ellipticity, and there- 
fore lead to stable material models. The standard fractional-order viscoelastic (fov) 
solid is used to arrive at the overall elastic/ viscoelastic structure of the model, while 
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the elastic potential — via the k-bkz hypothesis — is used to arrive at the tensorial 
structure of the model. Candidate sets of functions are proposed for the elastic 
and viscoelastic material functions present in the model, including a regularized 
fractional derivative that was determined to be the best. The Akaike information 
criterion (aic) is advocated for performing multi-model inference, enabling an ob- 
jective selection of the best material function from within a candidate set. 

Key words nonlinear elasticity — viscoelasticity — finite strain — fractional calcu- 
lus — Mittag-Lefller function — soft tissue — multi-model inference 


i Introduction 

The human heel is comprised of skin, a fat pad, the origin of the plantar aponeu- 
rosis tendon and the calcaneal bone. Collectively, the soft tissues therein constitute 
the heel pad. The heel pad is our body’s natural shock absorber, dissipating im- 
pulses introduced into the body during normal activity, and thereby attenuating 
the forces that are transmitted up through the body’s skeletal structure (Cavanagh 
et al. 1984). 

A certain amount of loading is needed to maintain healthy bone. Bed-rest stud- 
ies have demonstrated that prolonged inactivity depletes bone density (Shackelford 
et al. 2004). 

NASA has a need to understand how much force is being transferred into the 
load-bearing bones of the body during exercise so that effective countermeasure 
protocols can be developed to help avert bone loss in astronauts during long space 
missions (Lang et al. 2004). Current devices that measure in-shoe forces beneath 
the heel have recorded forces that exceed twice body weight when an astronaut ran 
on a treadmill on Earth; whereas, when running on an identical treadmill located 
within the International Space Station, using the same in-shoe transducers and 
a harness attached at the waist pushing the astronaut against the treadmill, maxi- 
mum forces of about one and one-halftimes body weight were recorded (Cavanagh 
et al. In press). To determine how much of this force is actually being transmitted 
to bone will require numerical analysis. To be able to run such an analysis will 
require material models for the soft-tissue constituents of the heel pad. Here we 
develop a viscoelastic material model for the human calcaneal fat pad. 

Although the fractional calculus 1 has enjoyed wide application in synthetic 
polymer rheology (see Podlubny 1999, pp. 268—277, for a brief literature review), 
it has attracted limited attention in the field of biomechanics: Suki et al. (1994) 
found the pressure/volume response of a whole lung to be aptly characterized by 
a Newtonian fov material model with a fractional order of evolution of o.i; 2 his 

1 Calculus is the study of properties of functions in one or more variables, using derivatives 
and integrals. Fractional calculus extends the classic study of integer-order derivatives and 
integrals to include derivatives and integrals of non-integer order, e.g., d n y /dx n . 

2 For the Newtonian fov model, a fractional order of evolution equaling 1 is the limiting 
case of a viscous Newtonian fluid, while a fractional order of evolution equaling o is the 
limiting case of an elastic Hookean solid. 



A K-BKZ Formulation for Soft-Tissue Viscoelasticity 


3 


colleagues, Yuan et al. (1997, 2000), studied lung tissue and found its fractional or- 
der of evolution to be about the same, viz., 0.075; while Chen et al. (2004) applied 
the same model to agarose gels used for culturing tissues, especially cartilage, and 
found its value to be about 0.03. These are all values close to that of ideal elas- 
ticity, where the order of evolution is o. In a study of charge dynamics in protein 
molecules, Glockle and Nonnenmacher (1995) derived a kinetic equation in the 
form of a fractional-order integral equation (i.e., a Volterra equation of the second 
kind with an Abel power-law kernel) and found the charge relaxation in myoglobin 
to be accurately described by a formula where the fractional order of evolution was 
set at -0.4. 3 In papers by Carew et al. (2003) and Doehring et al. (2004), the response 
of aortic heart valves to ID experiments has been shown to be well represented by 
a quasi-linear 4 Kelvin-Zener fov solid with the fractional order of evolution being 
about 0.25. 

In this paper, we present a k-bkz viscoelastic model tailored to the response 
of the human calcaneal fat pad loaded in compression. The paper begins with a 
presentation of the kinematic fields needed to construct such a theory. A novel class 
of nonlinear elastic solids is then presented that has great potential in the model- 
ing of soft tissues. We use the aic information theoretic (Burnham and Anderson 
2002) to select a ‘best’ model for the heel pad based on the compression data of 
Miller-Young et al. (2002) — a power-law model was found to be best. A review of 
ID fov restricted to infinitesimal strains follows, which is then transformed into 
an equivalent Boltzmann integral equation with a memory-function kernel. This 
formulation is useful in that it permits an extension into 3 D by applying the k-bkz 
hypothesis to construct a finite-strain theory. This hypothesis employs the elastic 
strain-energy function (previously selected) to establish the tensorial structure of 
the viscoelastic model. Like our elastic material class, there are many candidate 
models that belong to our viscoelastic material class. In addition to the fov kernel, 
four other kernel functions are considered as candidate models. The aic informa- 
tion theoretic was then used once again to select a ‘best’ model, this time based on 
the stress-relaxation experiment of Miller- Young et al. (2002). A regularized frac- 
tional derivative (rfd), introduced in App. B, was found to be the best viscoelastic 
kernel. 

Appendices are provided that: A, describe how to use the aic information theo- 
retic; B, list the candidate viscoelastic kernels; C, derive our viscoelastic constitutive 
formula using converted tensor fields, and establish how such formulae map into 
Cartesian space; and D, list series expansions that allow finite-strain constitutive 
formulae to be recast as infinitesimal-strain constitutive formulae. 


3 Fractional-order integration and differentiation can be defined as a single operator that 
is continuous over the order parameter; hence, the term differ-integration (Oldham and 
Spanier 1974). The accepted notation employs a minus sign (e.g., -0.4) when designating the 
order of integration, and a plus sign (viz. , 0.4) when designating the order of differentiation. 

4 A viscoelastic model is said to be ‘quasi-linear’ if: the linear strain (or forcing function) 
of classic (linear) viscoelasticity is replaced by a nonlinear strain measure, the kernel (or 
memory) function depends solely on time (i.e., strain-time separability applies), and only a 
first-order integral over time appears in the model. K-bkz models are quasi-linear. 
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2 Kinematics 

Consider a rectangular Cartesian coordinate system with orthonormal base vectors 
e\, e 2 and C3. We focus our attention on a mass point originally located by the set 
of coordinates X = (Jf| , X2, A3) assigned at an arbitrary reference time to in 
this coordinate frame. At current time t, this mass element is located by a different 
set of coordinates X = (x\ , X2, X3) in the same coordinate frame, while at some 
intermediate time — say s, to < s < t — it had coordinates x = (Xi ■ X2< Xi)- 

It is supposed that the motion of this mass point through space can be described 
by a one-parameter family (in time) of locations considered to be continuous and 
sufficiently differentiable to allow the following deformation gradients to be defined 

3 Xi - „ dxi ~ 3 Xi 

1-U(s.n:= w . F ij{ to , S ):=JL, (I) 

where indices i and j have values 1 , 2 , 3 . Here these formulas have been written in 
component form; in tensor notation they are written as 

F = Fjj ei <g> ej, F = Fy e,- <g> ej, F = Fy ?; <g> ej, (2) 

where ® is the vector outer product. These fields satisfy the identity F = F F , or 
equivalently, Fij = FjfcFkj where the repeated k index is summed over in the 
usual manner. The ability to invert these fields, guaranteed by the conservation of 
mass, ensures that a given particle cannot occupy two locations at the same instant 
in time, and that two discrete locations do not associate with a single particle at any 
given moment in time. 

Deformation fields are two-state fields that can be scalar, vector or tensor val- 
ued. Hereafter, arguments denoting the state dependence of these fields are omitted 
for brevity, at least for the most part. Instead, as in Eq. (2), plain-symboled deforma- 
tion fields are considered to have a state dependence of (to, t); hatted deformation 
fields are considered to have a state dependence of (s,t); and tilded deformation 
fields are considered to have a state dependence of (to, s). 

Affiliated with the above deformation gradients are the left- and right-defor- 
mation tensors defined by 


B:=FF t and C:=F t F, (3) 

respectively, where tT ’ implies transpose (viz., B t j = F^Fj^ and Cjj = FkiF^j). 
By B we mean ff t , etc. The left-deformation tensor B of Finger (1894) typically 
appears in Eulerian constructions, while the right-deformation tensor C of Green 
(1841) typically appears in Lagrangian constructions. 

For model implementation into numerical codes, like finite elements, it is of- 
ten useful to split the deformation variables into hydrostatic and deviatoric parts. 
Following Flory (1961), we assign 


J := detF, F := /“ 1 / 3 F, C := F T F, B := FF T , 


( 4 ) 
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so that det F = 1 , and therefore, det C = det B = 1 , where det(-) denotes the 
determinant. Likewise, one can define 

J := detF, F := /“ 1/3 F, C := F T F, B := FF T , (5) 

and _ _ _ _ _ 

J := detF, F := /“ 1/3 F, C := F T F, B := FF T , (6) 

so that det F = det F = 1 . A bar over a tensorial deformation field implies that it 
is isochoric (preserves volume), while a hat or a tilde on top of that designates the 
state in which it is isochoric. 


3 Elasticity 

Before one can construct a viscoelastic model for soft tissues, it is necessary to 
quantify the highly nonlinear elastic behavior that dominates soft-tissue response. 
Specifically, what we are after is the mathematical form of the elastic strain-energy 
function for the material of interest; herein, the human heel pad. Different tissues 
are typically governed by different strain-energy functions. 

The strain-energy density per unit mass, when written for the Lagrangian 
frame, is given by 

dW= -^tr(SdC), ( 7 ) 

where tr(-) is the trace operator, while dW(X', to, t, d?) represents the work done 
over a time increment d? on a material element with mass density Q = q(x\ t ), 
where Qo = Q{X; to). Work is caused by an imposed displacement acting on the 
mass element, manifested here as the strain increment ydC ( X ; to , t, d/). The ma- 
terial responds to this displacement through the creation of forces, thereby produc- 
ing a state of stress S(X \to,t). It follows that = det F from the conservation of 
mass. 

Elastic states occur at minima in the strain-energy function W. Therefore, an 
incompressible elastic solid is defined in the Eulerian frame by the constitutive law 

dW(C) r 

T + p I = 2 pF — — — F t , detF = 1, (8) 

0 C 

where the isochoric constraint for incompressibility, det F = 1 , is introduced into 
the formula through a Lagrange multiplier p multiplying the identity tensor I . Soft 
tissues are comprised primarily of water, whose bulk modulus is 2.2 GPa. Most of 
these tissues have shear moduli that typically range between a kPa and 10’s of 
MPa. Consequently, the ratio of their bulk to shear moduli usually lies between 
100 and 100,000, depending on the tissue, making incompressibility a reasonable 
assumption to impose for these tissues. An application of the push-forward operator 
(Holzapfel 2000, pp. 82—84) transforms the Lagrangian form of this law, which one 
obtains by minimizing Eq. (7), into the Eulerian expression presented in Eq. (8). 
The second Piola-Kirchhoff stress S maps into the Cauchy stress T (jc; t) according 
to the well-known formula T = ^-FSF T . 
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There are two non-trivial invariants needed to describe an isotropic, incom- 
pressible, elastic solid (Rivlin 1948); they are: 

/ = tr C = tr B & II = i((trC) 2 -tr(C 2 )) =trC _1 = trB -1 . (9) 

Consequently, the tensorial dependence of W(C ) can be replaced by a scalar one of 
W(I , H). The strain-energy function W will also depend on a number of material 
constants. The third invariant is a trivial argument, because IH = detC = 1 
from the incompressibility assumption. We will reintroduce the third invariant in 
§5 . An application of the C ayley-Hamilton theorem proves the identity j ( (tr C ) 2 — 
tr(C 2 )) = trC~* provided that det C = 1 , which follows from incompressibility. 

We are not free to assign any functional form of our choosing to the strain- 
energy function. Material stability needs to be taken into consideration. This is 
especially important so as not to inadvertently incite numeric instability into soft- 
ware codes like finite elements. A critical hypothesis needed for well-posedness of 
such initial value problems is that the equation of motion be hyperbolic or, equiv- 
alently, that the corresponding quasi-static constitutive equation be elliptic. In this 
regard, Renardy (1985) proved a lemma that establishes a sufficient condition for 
strong ellipticity in a k-bkz fluid that has relevance in elasticity as-well-as in k-bkz 
viscoelasticity whenever both invariants / and II are active. 


Lemma 1 A sufficient condition for strong ellipticity in an incompressible, isotropic, elastic solid 
is that its strain-energy function be strictly: monotone in I and II , and convex in \fl and \fH . 


The monotonicity part of Renardy’s lemma mandates that the gradients of the 
strain-energy function (taken with respect to / and II) must be positive; that is, 


dW 

~JT 


> 0 


and 


dw 

JF 


> 0 . 


(10) 


The convexity part of his lemma mandates that the Hessian of the strain-energy 
function (taken with respect to VT and cfll) must be positive definite, which via 
Sylvester’s theorem (a.k.a. the criterion ofHurwitz) requires that 


d 2 W 3 2 W 3 2 W ( 3 2 W V 

> 0 and — > I — — — 1 . (11) 

(3V7) (dVl) (dVF) \dVTdVTlI 

An important example of a strain-energy function that satisfies this lemma is the 
sum W oc I + II . Simple counterexamples that do not satisfy this lemma include: 
the product W oc I II , which is monotonic but not convex; the ratio W oc / / II , 
which is neither mono tonic nor convex; and the product sum W oc I 2 II + I II 2 , 
which is monotonic but not convex in a large region surrounding the stress-free 
state where B = I . Many other counterexamples can also be constructed. 

Rivlin and Saunders (1951) introduced a phenomenological class of incompress- 
ible elastic solids whose strain-energy function is given by 


00 

2 qW M r= C aP (I-3) a (n-3f, C oo = 0, (12) 

a=0J=0 
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where the C a p are material constants to be obtained through parameter estima- 
tion procedures. The literature often refers to this as the Mooney-Rivlin solid. The 
neo-Hookean solid (derived from statistical molecular physics) has a non-zero coef- 
ficient for C 10, while the phenomenological Mooney solid has non-zero coefficients 
for C 10 and Cot (Treloar 1975, pp. 212— 213). Both are popular models in the rubber 
elasticity literature. To satisfy Lemma 1, the coefficients in the neo-Hookean and 
Mooney solids must all be positive valued. Mixed coefficients C a p, where a > 0 
and P > 0, associate with terms that do not satisfy Lemma 1, in general. 

We propose an alternative phenomenological class of incompressible elastic 
solids that satisfies Lemma 1, and which finds application in biological tissues. El- 
ements of this material class have gradients of strain energy that produce formulae 
for stress T that are proportional to strain fields (e.g., ^(B — B^ 1 ) in the simplest 
case). It is the tensorial strain measure that becomes 0 whenever B = I in our 
material class. In contrast, the Mooney-Rivlin class of elastic solids in Eq. (12) pro- 
duces formulae for T that are proportional to the deformation fields B and — B -1 , 
and as such, they rely on the Lagrange multiplier p in Eq. (8) to absorb all resid- 
ual strain occurring at B = I. This is permissible, mathematically; it is just not 
desirable, physically, in our opinion. 

The Eulerian strain tensor j (B-B ') has a symmetry in its dependence upon 
state. It is also a second-order accurate approximation to the true strain measure 
of Hencky, viz., j In B (Freed 2004). Contrast this with the un-symmetric Eulerian 
strain measures ofSignorini (1930) 4 (B — I) and Almansi (1911) ^-(I — B _1 ) that are 
prevalent in the literature, both of which are first-order accurate approximations 
to Hencky strain. Einstein et al. (In press) and Freed et al. (In press) have found 
j(B — B" 1 ) to be a good strain measure for quantifying the response of the ground 
substance matrix in soft-tissue mechanics, where they employed this Eulerian strain 
measure in its Lagrangian representation which reads j(I — C^ 2 ). 

Whenever an invariant appears by itself in a strain-energy function, a defor- 
mation tensor ensues. Whenever the invariant sum / + II appears, a strain tensor 
is produced. Our material class is a generalization of this invariant sum, keeping in 
mind the constraints of Lemma 1 . Let us consider a class of incompressible elastic 
materials whose strain-energy function is given by 

2qW = - f(pu 3) + f(p 2 \ II) - /O 2 ; 3)), ( 13 ) 

where (> 0 ) is the elastic shear modulus with units of stress, and p\ is a vector 
of parameters, which may have different values when associated with / and II . 
Function f is any dimensionless function that satisfies the following constraints: 

physics: f(pi', /) > 0 & f(pi\ H) > 0, 

f'ipu 3)= 1 & f\ P2 \ 3) = L 

monotonicity : f \p\ \ I) > 0 & f'{ P2 -n)>Q, (i 4 ) 

convexity: f'(Pi’I) + 2If"(j>i',I) > 0 & 

f\pi, n) + 2 nf"( P2 ; n) > o, 

where f'(pc) := d/(A')/d.v and f"(pc) := d 2 f(x)/dx 2 . The /(/>,•; 3 ) present in 
Eq. (13) are constants introduced to normalize the strain energy so that W > 0. 
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The definition for strain energy given in Eq. (13), in conjunction with the defin- 
ing law for incompressible elasticity stated in Eq. (8), lead to an elastic constitutive 
equation in the Eulerian frame of 

T detF = 1 , (15) 

that when pulled back into the Lagrangian frame becomes 

S + pC- 1 = f(p2\n)C- 2 ), detF = 1 , (16) 

recalling that Q = Qo because detF = 1 . The V4 is introduced here so that 
whenever f = 1 the above equation reduces to T + p I = /x^B—B -1 ), wherein 
fl now corresponds with the classic definition of Lame’s elastic shear modulus in 
the domain of infinitesimal strains. This physical interpretation of fl applies to all 
models in our material class, because of the second line of constraints in Eq. (14). 
These same constraints also ensure that the right-hand side of Eq. (15) becomes 0 
whenever B = I, i.e., the right-hand side is a strain measure independent of p. 

Choosing a functional form for f that is in accordance with the constraints put 
forth in Eq. (14) will lead to an admissible constitutive equation for the modeling of 
elastic solids. 


3.1 Human Heel Pad 

Few tissues in the human body are isotropic. The calcaneal fat pad in our feet 
has been demonstrated via experiment to be isotropic (Miller- Young 2003). This 
makes the heel pad an ideal tissue to work with for the purpose of deciphering 
mathematical structure, and to assess capability of the aic information theoretic in 
model selection through multi-model inference — a technology reviewed in App. A. 
Extending the structure of our model to anisotropic tissues is a topic for future 
work. 

There are a variety of functional forms for f that one could investigate which 
will satisfy the constraints laid down in Eq. (14). We shall consider four models 
constructed from the following two basic mathematical functions: 

f(x)= T ^ T x n+ \ n> 0, xe {1/3, E/3}, 

f(x ) = U nx , n > 0 , x€{I-3,E- 3 }. 

Parameter fl is common to all models, while the parameter vectors pi are equal 
(i.e., pi = pi = {«}) in two of the models, and distinct (viz., p \ = {n 1} and 
Pi = {/I2}) hi the other two models. The power law has a long history in tissue 
mechanics, dating back to Mitton (1945). More prominent in the biomechanics 
literature of today is the exponential law advocated by Fung (1967). 

We proceed by acquiring maximum log-likelihood estimates for the unknown 
material parameters in each of the four candidate models, along with a suite of 
statistical parameters: the objective function <?, the coefficient of variation in the 
data o', the aic information theoretic fi Alc and the aic difference A;, as defined in 
App. A. These values have been tabulated in Table 1. Quasi-static and dynamic 
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fix) 

Mqs (kPa) 

M'dyn (kPa) 

n 1 

n 2 

0 

O 

M AIC 

A, 

x n 

0.691 

2.214 

2-59 

n.2 = w 1 

5.5662 

0-2395 

1 75 -° 

0 

x n 

p 

00 

C7I 

2.177 

5.08 

1.29 

5 - 54 I 4 

0.2390 

776.7 

!-7 

e nx 

O.7OO 

2.274 

O.708 

«2 = n 1 

5.6306 

0.2409 

176.1 

1. 1 

e nx 

O.692 

2.227 

I.4O 

0.405 

5-5722 

0-2397 

T 77-3 

2-3 


Table i Optimized parameters /r qs , /tdyn, 7*1 and «2 for the quasi-static and dynamic 
elastic responses of the human calcaneal fat pad in unconfined compression, cf. Figs, (i & 2). 



perimental mean and standard deviation data (obtained from 10 feet) are from Miller- Young 
(2003). 


experiments were fit simultaneously — see Figs. (1 & 2). It was postulated and ver- 
ified that only the shear modulus /x exhibits a rate dependence, whereas n is rate 
insensitive in a statistical sense, which is why there are two shear moduli reported 
in Table 1; they are the shear moduli belonging to these two experiments, and are 
not to be confused with the viscoelastic rubbery /Too and glassy /To shear moduli 
introduced in the next section. 

By employing the methodology from information theory presented in App. A, 
an examination of the data presented in Table 1 allows one to conclude that the 
power-law and exponential models are both ‘good’ candidates for the modeling of 
unconfined compression in the human heel pad, with the power-law being only 
slightly better in this instance. The power law has an additional practical advan- 
tage over the exponential in that it is more efficient and robust in a finite element 
setting. For both function types, the models with n \ = 112 were found to be supe- 
rior to their affiliated models where n \ 7^ 77 2. The additional parameter present 
in the models where n 1 7^ 7? 2 brought no added value to these models from the 
perspective of information theory, allowing the simpler models where n\ = no to 
be selected. Herein lies the true worth of the aic information theoretic: models 
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Figure 2 Dynamic stress/stretch response to 50% deformation at X = —35 s 1 . Exper- 
imental mean and standard deviation data (obtained from 7 feet) are from Miller- Young 
(2003). 


with differing numbers of material parameters can be assessed objectively to deter- 
mine which is best. Aic provides a metric for the Kullback-Leibler (kl) information 
space by which distances can be measured between various mathematical models 
(via the aic differences A{) whose parameters are all fit against a common data set. 

Fits of the power-law model (using Eqs. 15 & 17 with n \ = 722) to the quasi-static 
and dynamic, experimental data sets of Miller- Young et al. (2002) can be found 
in Figs. 1 & 2. Although not perfect, these fits are within experimental variation 
up to about 45% compression (A. = 0 . 65 ). There is too much curvature in the 
model to correlate the quasi-static data with exacting precision. In contrast, there 
is not enough curvature in the model to accurately correlate the dynamic data. The 
model should therefore provide reasonable approximations of reality over a wide 
range in dynamic input. Placing 95% confidence intervals around the parameters 
of this fit (see Eq. A4) puts n G [ 2 . 30 , 2 . 85 ], while /r qs G [ 0 . 48 , 0 . 90 ] kPa and ix&yn G 
[ 1 . 55 , 2 . 89 ] kPa. These are maximum-likelihood confidence intervals, which need 
not be symmetric about their optimum values, as is the case with least-squares 
confidence intervals. 


4 ID FOV 

In landmark papers by Caputo and Mainardi (1971a, b), the authors analytically 
continued the standard viscoelastic solid (Zener 1948, pg. 43) 


[l + r D]ct(?) = T'oofl +pD]c(f), ct 0 + = E oo (p/r)e 0 +, 


(18) 
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by replacing its derivatives in time D f(t) := df(t)/dt with the Caputo (1967) frac- 
tional derivative of order a in time (cf. Podlubny 1999, pp. 78— 81) 5 

1 C‘ D f(s) 

:= TTTTi) / 0 (F^sr ds> C9) 

thereby producing the constitutive equation 

[l + t“D“]o-(» = ^oofl + /o“D“]e( 0 , or 0 + = E oo (p/r)°‘e 0 +, (20) 

that we call the standard fov solid 6 , which becomes the Kelvin-Zener viscoelastic 
solid listed in Eq. (18) whenever a = 1 . Variables er and e represent engineering 
stress and strain, respectively, with cr 0 + := cr(? 0 +) and e 0 + := e(t 0 + ) specifying 
their initial conditions at time t 0 + (= to + £, where £ is a small positive number). 
This ID model has four material constants: E 0 a (> 0 ) denotes the rubbery elastic 
modulus, a (0 < a < 1) is the fractional order of evolution, r (> 0) represents the 
characteristic relaxation time, and p (> r) is the characteristic retardation time, 
with Eq := (p/x) a Eoa (> E 0 o) establishing the glassy elastic modulus. 

Bagley and Torvik (1986) have shown that the fractional orders of differenti- 
ation for stress and strain must be equal, as originally proposed by Caputo and 
Mainardi (1971a), in order to ensure that this constitutive relationship is compati- 
ble with the second law of thermodynamics; specifically, they must be equal so as 
to guarantee non-negative dissipation during cyclic loadings. 

Having unequal fractional derivatives on the two sides of the equation also in- 
troduces unbalanced shocks into the solution arising from the initial conditions, 
and is most evident in the Laplace domain (Bagley and Calico 1991). Shocks need 
to be in balance in order for the predicted stress waves to travel with finite veloc- 
ity, in accordance with physical observations. Not all fov models possess balanced 
shocks; for example, the popular Voigt fov solid ar(t) = 1 + p“D“]€(f) intro- 

duced by Caputo (1967) has an unbalanced shock, and as such, predicts that stress 
waves will travel with infinite velocity. Material models whose stress waves are pre- 
dicted to travel with infinite velocity are a known source of numeric instability in 
finite element codes that account for inertial effects (Belytschko et al. 2000, pg. 
3 H)- 


5 In App. B we introduce a kernel that regularizes the fractional derivative so that Eq. (ig), 
for 0 < a < 1 , takes on the form 


D %f(f) 


m-a) J t0 


p m 

(t + 8 - s) a 


dj. 


8 > 0, 8 /(t — t 0 ) « 1, 


where 8 is a small positive number (relative to t ) that effectively takes the singularity at the 
upper limit of integration in the Caputo derivative (19) and moves it a minute distance 8 
outside the interval of integration. See §4 in App. B for more details. 

6 The standard fov fluid is defined by 


[1 + r“D"]<r(0 = rj a D“e(t), a 0 + = (77/ r)“e 0 +, 


where tj (> 0) is the viscosity, r (> 0) is its characteristic relaxation time, and a (0 < a < 1) 
is the fractional order of evolution. 
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A robust and efficient numerical algorithm capable of solving the fractional- 
order differential equation presented in Eq. (20) has been published by Diethelm 
et al. (2002, 2004, 2005). 


4.1 Relaxation-Function Formulation 

An analytic solution to the standard fov solid (Eq. 20) was obtained by Caputo 
and Mainardi (1971a) through an application of the method of Laplace transforms. 
They are able to apply this technique to Eq. (20) because it is a linear differential 
equation, albeit of fractional order. The solution they arrived at is a special case 
of Boltzmann (1874) viscoelasticity, commonly written as (cf. Christensen 1971, pp. 
3 - 9 ) 

o(t) = Q{t) e 0 + + f Git - s ) d.s. (21) 

Scaling the relaxation function G(t) so that it reads as 

Git) = £00 + (T 0 - Too) Git), (22) 

with the normalized relaxation function G)t) constrained so that Go := G( 0 ) = 1 
and Goo •= G(oo) = 0 , 7 allows one to rewrite Eq. (21) as 

cr it) = E^it) + (T 0 - Too) ^G(t) e 0 + + J* Git-s) d.v j , (23) 

which will describe a viscoelastic solid if To > Too > 0, and a viscoelastic fluid if 
Eo > Too = 0 with the lower limit of integration then set at to = —00. 

For the standard fov material models, the relaxation function has the special 


form 

Git) = E aA (-{t/x) a ), 

( 24 ) 

with 

00 z k 



E a,p(z)- 'E r(J} + ak y Z€C, 

( 25 ) 


defining the Mittag-Leffler function (cf. Podlubny 1999, pp. 16—37), wherein R de- 
notes the real line, R+ the positive real line, and C the complex plane. Equation 
(24) satisfies the constraints Go = 1 and Goo = 0 . The Mittag-Leffler function first 
appeared as a relaxation function (Eq. 24) in a paper written by Gross (1947), where 
it was introduced in an attempt to remedy inconsistencies present in the power-law 
creep function. Gross did not connect the Mittag-Leffler relaxation kernel with the 
fractional calculus. That took place later in the papers of Caputo and Mainardi 
( I 97 la 5 b). 

Relaxation functions described in terms of the Mittag-Leffler function, as occur 
in phenomenological fov models, arise naturally from the statistical mechanics of 


7 A relaxation function with G 0 = 00 is indicative of a material that propagates an 
impulse with infinite speed (like the Voigt fov solid), which is not physical. 
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random walks made with steps taken at random intervals (Douglas 2000). Other 
examples of relaxation functions that have been used in soft-tissue mechanics are 
listed in App. B. 

A straightforward summation of the series presented in Eq. (25) has regions 
of numeric instability whenever computations are done on machines with finite 
precision; therefore, equivalent expressions need to be employed in various sub- 
domains of the function space to ensure an accurate computation of E a jj (z) for 
all admissible values of a, and z. Gorenflo et al. (2002) have published one such 
algorithm, which is more readily available in the paper of Diethelm et al. (2005). 

When the fractional order a is very small, but still greater than zero, the relax- 
ation kernel G(t ) = E a y(—(t /r)“), plotted in Fig. 3 for r = 1 , drops immediately 
from a value of G = 1 at / = 0 to a value of G ss */2 at t = 0 + , from which it 
asymptotes algebraically, and very slowly, to zero as t goes toward infinity. Strictly 
speaking, the Mittag-Leffler function is not defined at a = 0 . This would be the 
elastic boundary where, if the Mittag-Leffler function were defined, it would ex- 
hibit a discontinuous jump so that at a = 0 the response would be G = '/ 2 at 
t = 0 and G = 1 for all t > 0 (i.e., the Heaviside unit-step function), as there is 
no relaxation in the elastic limit. At the other boundary where a = 1 , relaxation 
is smooth and asymptotes exponentially to zero as t moves toward infinity because 
E\,\(-t/ r) = e~ f/T . 

For all values of a that lie within the interval ( 0 , 1 ], the Mittag-Leffler function 
provides a smooth monotone transition from a perfect elastic to a classic viscoelastic 
response. The function is not monotone whenever a > 1 . 


4.2 Memory-Function Formulation 

Infinitesimal strain is actually a two-state field that we can write as e(to,t), where 
time to denotes some reference state, and time t represents the deformed state. 
After an integration by parts, Boltzmann’s viscoelastic model (23) becomes 

ar(t) = E 0 e(t 0 , t) - ( E 0 - E 0 c ) f M(t-s) e(t 0 , s) d.v, (26) 

J to 

wherein M(t—s) := dG(t — s)/ds defines the memory function used by Lodge 
(1956) in a model that he latter called the rubberlike liquid (Lodge 1964, pp. 101— 
104). The notion of a memory kernel was re-introduced into the biomechanics 
literature by Zhu et al. (1991). The terminology ‘memory function’ is a synonym 
for the ‘rate-of-relaxation function’. 

By utilizing additivity of infinitesimal strains (i.e., e(?o, t) = e(?o, s) + e(s, t) 
for all s e [/o, ?]), the above formula can be recast as 

ct( 0 = ( E 00 -\-(Eo-E 00 )G(t))e(to,t) + (Eo-E 00 ) f M(t-s) e(s,t)ds, (27) 

Jt 0 

where now the reference state is s in the strain variable that lies under the integral 
sign. This is consistent with the physical notion that interval [v , t\ constitutes that 
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part of the deformation history which the material recollects, while the preceding 
interval [?o, s) represents that part of the history which the material has forgotten. 
In this regard, Eq. (27) is preferred over Eq. (26) because it better reflects the un- 
derlying physics, even though Eq. (26) has a simpler structure and the two formulae 
are equivalent. Equation (27) is also in a form amenable to the k-bkz hypothesis. 

In the testing of soft tissues, it is customary to precondition the sample prior 
to executing a test (Fung 1993, pp. 260—262). This has the effect of rendering 
down the non-integral viscoelastic contribution present in Eq. (27) to its quasi-static 
constituent, viz., E^eito, !)> an d as such, 


o(t) = Eooe(to,t) + (E 0 - Eoo) [ M(t-s)e(s,t)ds (28) 

J t 0 


becomes the governing constitutive expression for preconditioned specimens. 

The constitutive formulae in Eqs. (26—28) are less restrictive than the constitu- 
tive formula in Eq. (23). Equations (26-28) require strain to be a C 1 function (i.e., 
continuous and differentiable over the deformation history); whereas, Eq. (23) re- 
quires strain to be a C 2 function (viz., continuous and twice differentiable over the 
deformation history). 

Memory fades if 0 < < M(t\) for all ?2 > t\ > to, which is actu- 

ally a thermodynamic requirement (Coleman and Mizel 1968). This constraint is 
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Figure 4 A 3D plot of the fov memory function (i.e., Eq. 29 ), M(t) = 
— E a ,o(— (f/r) a )/f, with r = 1 where t € [0.001, 3], 


satisfied by the memory function for the standard fov solid 


M(t) = 


t 


(29) 


where, notably, E a fi(x) appears in the memory function, while E a \{x) appears 
in the relaxation function. The derivative dE a jj(x)/dx, which is required because 
of M(t — s) = dG(t — s)/ds, can be found in Podlubny (1999, pg. 22), for exam- 
ple. This kernel also appeared in the paper of Gross (1947), but it was written as 
—dE a (—(t/r) a )/dt, where E a (x) = E a ^\(x) is the one-parameter Mittag-Leffler 
function. Gross did not make use of the two-parameter Mittag-Leffler function 

E(i,fl(x)- 

When the fractional order a is very small, but still greater than zero, the mem- 
ory kernel M(t) = — E a- o(— (t /z) a )/t, plotted in Fig. 4 for r = 1 , has a response 
M that behaves like an impulse function, indicating that the material has a per- 
fect knowledge of the current state. In contrast, it has virtually no recollection of 
even the most recent of past states. As a approaches unity, the memory function 
continues to maintain its perfect knowledge of the current state (i.e., M is infinite 
at t = 0, except at a = 1, however the strength of this singularity diminishes as 
a — > 1 ). To this complete remembrance of the current state, the function then adds 
an increasing recollection of past events with increasing a; albeit, this is a memory 
that fades away rapidly with the passage of time. 

The fact that the memory function in Eq. (29) possesses a weak singularity at 
Mq := M{ 0 ) for all values of a € ( 0 , 1 ) needs to be taken into consideration when 
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Figure 5 A 3 D plot of M(t)/e 1 = —E ol 0 (—{t/r) a )/(te f / r ) with r = 1 , where t e 
[ 0 . 001 , 5 ], demonstrating that the fov memory function in Eq. (29) algebraically asymptotes 
to zero as t — > 00 for all a 6 (0, 1). 


selecting a numerical method for solving the convolution integral that appears in 
Eqs. (26—28). We have constructed a unique memory-management scheme, and 
have coupled this with a midpoint quadrature rule with a Laplace end correction 
to produce an efficient (i.e., 0 (N log N)) and accurate (viz., 0 (h 5 )) numerical al- 
gorithm that solves convolution integrals like those listed in Eqs. (27 & 28), wherein 
the forcing function depends on both states s and t, instead of on just state t as is 
usually the case in convolution integrals (Diethelm and Freed in review). 

A visual inspection of Figs. 3 & 4 indicates that there should be a significant nu- 
merical advantage when employing the memory function defined in Eq. (29) over 
its corresponding relaxation function given in Eq. (24) for the kernel of viscoelastic 
convolution, provided that the singularity at the upper limit of integration can be 
effectively and efficiently handled. Further inspection of Fig. 4 may mislead one to 
draw a false conclusion that the standard fov memory function fades faster than 
the exponential, which is refuted in Fig. 5. For a € ( 0 , 1 ) and an argument t/x 
that is less than about three, the fov memory function does indeed collapse faster 
than exponential decay, except in a neighborhood around the origin. Flowever, as 
the argument exceeds three, this trend begins to reverse and the fov memory func- 
tion starts to exhibit its true character of being algebraically asymptotic. Decay is 
exponential only when a = 1 . 


Remark 1 The integral equation given in Eq. (27), when employing the material 
functions defined in Eqs. (24 & 29), is an equivalent representation of the fractional- 
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order differential equation stated in Eq. (20): the standard fov solid. A practical 
advantage of expressing FOV as a Boltzmann integral equation with a 
Mittag-Leffler kernel, instead of as a fractional-order differential equa- 
tion, is that a working knowledge of the fractional calculus is not re- 
quired of the user in order for him/her to understand FOV to the extent 
that he/she can use it with confidence to solve problems of engineering 
interest. Of course , this analog only exists if the viscoelastic response is linear for 
the material of interest. 


5 Viscoelasticity 

We now extend this ID formulation, and in particular Eq. (27), into a 3 D theory 
that can be used to model soft isotropic tissues. To achieve this objective, we em- 
ploy the k-bkz hypothesis. 8 This hypothesis takes the potential structure for elastic- 
ity arising from thermostatics and analytically continues it into neighboring states 
of irreversibility where viscoelastic phenomena exist. The thermodynamic admis- 
sibility of this hypothesis is discussed in a separate paper by Bernstein et al. (1964). 

Because most tissues are predominantly elastic, with a secondary viscoelastic 
attribute, and because the elastic response in these tissues is highly nonlinear, we 
believe that the k-bkz hypothesis has an advantage over other approaches when it 
comes to developing viscoelastic models for soft tissues, the most notable alternative 
approach being that of internal state-variable theory (Coleman and Gurtin 1967). 

Our viscoelastic model is actually constructed in App. C using converted tensor 
fields. The resulting formulae are then mapped into Cartesian space-tensor formulae 
by using transformation rules that are also presented in App. C. The outcomes 
of these transformations are the objective material models presented below: one 
each for the Eulerian and Lagrangian frames of reference, and a third that applies 
whenever strains are infinitesimal, whose derivation relies on the series expansions 
given in App. D. 


5./ Lagrangian Formulation 

A viscoelastic material model has been derived in terms of converted tensor fields 
in App. C by employing the k-bkz hypothesis. Transforming this converted model 

8 Bernstein et al. (1963) state their hypothesis thusly: “For the Coleman-Noll fluid, the 
stress at time t depends upon the history of the relative deformation between the configu- 
ration at time t and all configurations at times prior to t . To this idea we add the following 
notions: (1) The effect of the configuration at time r < t on the stress at time t is equivalent 
to the effect of stored elastic energy with the configuration at time T as the preferred config- 
uration. The effect depends on t — r, the amount of time elapsed between time r and time 
t. (2) The stress at time t is the sum (integral) of all the contributions from all x </.... In 
effect, we are taking the concept of a strain energy function associated with the theory of 
finite elastic deformations, which is formulated in terms of a preferred configuration, and 
incorporating it in a fluid theory of the Coleman-Noll type by treating all past configurations 
as preferred configurations.” 
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into Cartesian space in the Lagrangian frame, using mappings that are also pro- 
vided for in this appendix, leads to a decomposition in the second Piola-Kirchhoff 
stress S of the form 

S + //tC _1 = £, tr(£C) = 0 , (30) 

where the hydrostatic pressure p = — y/ -1 tr(SC) and Lagrangian extra stress 
£ are given by separate constitutive formulae. The extra stress is deviatoric in the 
sense that tr(£C) = 0 . Here we do not impose a constraint for incompressibility, 
recalling that = J = det F = V det C = \flE from the conservation of mass. 
Hydrostatic pressure is taken to be governed by the constitutive formula (C7) 

p = - k W~ /-1 )’ ( 3 1 ) 

wherein k is the bulk modulus. The bulk response is not considered to be viscoelas- 
tic in vivo in soft tissues, at least in a rate controlling sense. 

The deviatoric response is taken to be governed by the constitutive equation 


£ = 2(7100 + (/to - /too) G(t))J 2/3 

x DEV[i(/'Q> i: /) I - f'(p 2 \ n ) C- 2 )] + 2 (/t 0 - /too)/" 2 / 3 


x j‘ M(t—s) DEV[±(/'Q> 1; /) C- 1 


-f(p 2 J)C- l C 



(32) 


where DEV[-] := (•) - \ tr [QC] C” 1 denotes the Lagrangian deviatoric operator, 
while /too and /to are the rubbery and glassy shear-moduli, respectively. Time to is 
associated with a stress-free equilibrium state. 

The relaxation G and memory M functions can be of whatever form one 
chooses (cf. App. B). They are not specified by the construction, only constrained 
in that M(t—s ) = 3 G(t—s)/ds, 0 < M(t 2 ) < M(t\) for all t 2 > t\ > to, Go = 1 
and Goo = 0 - Whenever G and M are given by Eqs. (24 & 29), respectively, Eq. 
(32) becomes a 3 D fov material model. 

The elastic function f is left unspecified by the general construction, too. It 
is, however, constrained by Eq. (14) so that the overall model satisfies the k-bkz 
stability criterion of Renardy (1985). 


§.2 Eulerian Formulation 

Mapping the convected model derived in App. C into Cartesian space in the Eu- 
lerian frame yields a decomposition in the Kirchhoff stress P(x; to, t) of the form 

p + jp i = n, tr n = 0 , (33) 

where the hydrostatic pressure p '.= — ytrT = — -j/^'trP and the Eulerian 
extra stress n are governed by separate constitutive formulas. The Kirchhoff stress 
P relates to the Cauchy stress T via the well-known identity P = ^-T so that 
P = FSF T . The Eulerian extra stress is deviatoric in an Eulerian sense in that 

tr n = 0. 
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The hydrostatic pressure p is still governed by Eq. (31). The Lagrangian extra 
stress Y of Eq. (32) pushes forward into the Eulerian frame (i.e., II = FHF T ) 
producing 


n = 2(mco + (mo - Moo) G( 0 ) d ev[±(f'(p i; /) B - f'(p 2 ; II) fi" 1 )] 


+ 2(/to — Moo) 


J' M(t—s) d ev[l (/'(/»! ; I) B - f(p 2 ; II) B _1 ) 


d.v, 


(34) 


wherein dev[-] := (•) — -j tr(-) I denotes the Eulerian deviatoric operator. 

The structure of the Eulerian model should be more intuitive than its La- 
grangian counterpart, at least for most people. 


5.3 Infinitesimal- Strain Formulation 

The theoretical constructions of this paper are geared for soft biological tissues. 
In the case of hard biological tissues, like bone, and many engineered materials, 
like plastics, infinitesimal strain analysis suffices. For these materials, the previous 
model simplifies to an expression for engineering stress a of the form 

a + p 1 = 1 , \tl= 0 , ( 35 ) 

where the hydrostatic pressure p = — ytr a is governed by 

p = —k tre, ( 36 ) 

and the deviatoric stress Z is described by 

E = 2(moo + (Mo — Moo) G(t)) dev[e] + 2(mo — Moo) [ M(t-s)&&\\e]ds, ( 37 ) 

Jt 0 

where e = e(to,t) denotes engineering strain, with e = e(s,t) = e — e given that 

e = e(to,s). 

Equations (35—37) were obtained by applying the approximation formulae listed 
in App. D to the above stated model in either its Lagrangian or Eulerian represen- 
tation, it does not matter which. 


5 .^ Human Heel Pad 

As an example of an isotropic viscoelastic tissue, we consider the human calcaneal 
fat pad; in particular, the experimental data of Miller- Young et al. (2002). The cri- 
terion of isotropy was experimentally verified, and an assumption of incompress- 
ibility has been imposed. Contrary to most soft-tissue testing, their experiments 
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were not preconditioned. The experiments were done in unconfined compression, 
and under these boundary conditions Eqs. (33 & 34) simplify to 

F IF 

T\\ = — = — — = 2(/ioo + (/x 0 - fi 0 o) G(?)) 

A A 0 

X /) (A 2 - A, -1 ) + f'(p 2 - II) (A - A- 2 )) + 20*0 - Moo) 

xf M(t—s) \ I) (A 2 — A -1 ) + f'{p2\ F)(X — A -2 )) d.v, (38) 

■I to 

where / = A 2 + 2 A -1 and II = 2 A + A -2 establish the two invariants under 
incompressible uniaxial loading conditions, with A being replaced by A in / and 
II , which are arguments in f. Here F is the applied force, Aq and A are the 
initial and current cross-sectional areas, and A = I/Iq and A = are the 

two stretches present, with lo, l s and l representing the initial, intermediate and 
current gage lengths of the specimen, respectively. 

Because the loading history was not recorded by Miller- Young (2003) for their 
stress-relaxation experiment, we were forced to impose an idealized loading history. 
This is not desirable (Doehring et al. 2004; Gimbel et al. 2004), but it is the best 
one can do in this case. A deformation rate of A = — 100 s _1 was assumed, which 
is in the vicinity of the uppermost capabilities of modern servo-hydraulic testing 
equipment. This rate was applied for 0.004 s to produce a final stretch of A = 0.6 
that was then held fixed for one minute. This loading history allows the integral in 
Eq. (38) to be decomposed into the sum of two integrals. The first integral is over 
the interval of loading 1 € [?o. ?i], and the second integral is over of the interval of 
relaxation t e For this particular experiment, to = 0 s, t\ = 0.004 s and 

?2 = 60 s. The advantage of breaking the integral into a sum of two integrals is that 
the second integral vanishes under the boundary conditions of stress relaxation, 
because A = 1 for all s € and therefore, the forcing function (viz., strain 

from s to t) is zero over the entire region [?i , fe]- All arguments t in the integrand 
of Eq. (38) remain t whenever t > t\. Only the upper limit of integration gets 
changed from t to t\ in the contributing integral. 

Following the method of approach used to select an elastic model, first a set of 
candidate viscoelastic kernels was chosen, and then the aic information theoretic 
of App. A was employed to down-select the better models at describing a speci- 
fied experimental data set; in this case, the stress-relaxation experiment of Miller- 
Young et al. (2002). The set of candidate models chosen for consideration includes: 
fov, gmm, kww, £LV and rfd, of which fov was detailed in §4 while the latter four 
kernels have been described in App. B. 

Because the loading data were not recorded, we assigned a value of 2.6 to the 
exponent n in the previously selected elastic power-law function f'(x) = x n of 
Eq. (17), where no distinction is made between p\ and pi , i.e., p\ = pi = {«}• 
This value is in agreement with our findings from fitting the elastic data, and with 
the observation that it is the shear modulus /X, not the strain exponent n, that ex- 
hibits rate dependence. In all models except gmm, this leaves four parameters to be 
obtained via parameter estimation techniques, whose values are listed in Table 2. 
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Model 

Moo (kPa) 

fi 0 (kPa) 

G 

c 2 

0 

0 

Maic 

A; 

FOV 

0.707 

4.04 

0.472 

0-389 

0.00143 

0.0114 

-50.1 

3- 1 

KWW 

0.921 

4.64 

0.263 

0.272 

0.00166 

0.0123 

-48.4 

4-8 

QLV 

0.965 

3-4 1 

0.0059 

3 1 - 8 

0.00229 

0.0144 

-44.9 

8 -3 

RFD 

0.711 

3-79 

0.366 

0.069 

0.00108 

0.0099 

-53-2 

0 


Table 2 Optimized shear moduli /too and )Iq, and viscoelastic kernel parameters denoted 
as c i and c 2 (see the body of the text for the mappings to their specific model parameters) 
for a stress relaxation of the human calcaneal fat pad, cf. Fig. (6). 


Parameters in common betwixt all five models include the rubbery [i 0 0 and 
glassy fj . o elastic shear moduli, and the elastic stretch exponent n = 2 . 6 . Except 
for gmm, each kernel has a relaxation/memory function pair with two material 
parameters that we denote as C\ and C2 in Table 2. These are not the notations 
that exist elsewhere in this paper, so here we establish their mappings and their 
units: For fov, c \ := a and C2 '■= r (s); for kww, c \ := and C2 := r (s); for qlv, 
Ci := ti (s) and C2 := t2 (s); and for rfd, ci := a and C2 := 8 (s). 

According to the selection criteria put forth in App. A, rfd is a ‘good’ model 
for inference, fov lies on the boundary between ‘good’ and ‘mediocre’, while both 
kww and qlv are ‘mediocre’ models in this regard for this material. Given this 
fact, rfd is the model of choice. The computational effort required to evaluate 
the rfd kernel is less than the computational effort required to evaluate any of 
the other kernels — an added bonus. The ability of rfd to correlate these data is 
demonstrated in Fig. 6. We reiterate that this selection process is based on the 
a priori assigned set of candidate models, and on the experimental data set chosen 
to fit. Different results are likely to follow given different materials, data sets and 
candidate models. 

This outcome of rfd being the ‘best’ model for inference came as somewhat 
of a surprise to us. Our personal bias going into this exercise would have been 
to select the fov kernel; this bias being based on many physically sound reasons. 
Biomechanicians would be apt to preselect qlv based on the biases of their back- 
grounds. The fact that qlv is not a good model for plantar soft tissue agrees with 
the recent findings of Ledoux et al. (2004). The capability of the rfd kernel, which 
is a generalization of the sifs kernel proposed by Johnson et al. (1996), and the ease 
by which it can be computed cannot be disputed. Other than around the origin, 
the rfd kernel is not all that different from the Abel kernel of the fractional deriva- 
tive present in the Voigt fov model, or the Mittag-Leffler kernel present in of the 
Kelvin-Zener fov model derived in §4, but it is a lot easier to work with. In effect, 
the rfd kernel slides the singularity at the upper limit of integration in the Voigt 
fov kernel so that it lies just outside the integral by a small distance of 8 . We coined 
the acronym rfd from the phrase regularized, fractional derivative , because it behaves 
like an Abel kernel whenever t 8 , but it does not propagate a shock wave with 
infinite velocity like the Voigt fov kernel does due to the regularization imposed on 
the rfd kernel, viz., Go = 1 for the rfd relaxation function. 
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Figure 6 Stress relaxation response at 40% deformation (A. = 0 . 6 ). Experimental mean 
and standard deviation data (obtained from 7 feet) are from Miller- Young (2003). The mean 
maximum true stress was -12.7 kPa. 


For the rfd model, placing 95% confidence intervals around the parameters 
puts: Hoc e [0.702,0.721] (kPa), /x 0 € [3.75,3.83] (kPa), a € [0.363,0.369] and 
8 € [0.066, 0.072] (s), with n fixed at 2.6 in accordance with our elastic findings. 
These confidence intervals are very tight when contrasted with those obtained for 
the elastic model. This is because of the high precision of fit attained with the relax- 
ation data, as contrasted with the more moderate fit achieved with the compression 
data. The above confidence interval for fj , lies within the confidence interval for 
/U.qs obtained in §3, implying consistency between the data sets. 

Conspicuously absent from the prior discussion is the Gmm model, which is 
the de facto standard of the viscoelastic literature at large. The number of Maxwell 
chains (i.e., mm elements) considered will affect the number of material parameters 
present in any given gmm model. It is not uncommon in the literature to find in- 
vestigators using upwards of 7 to 10 Maxwell elements in order to get a reasonable 
fit to a given set of experimental data. Nowhere, to our knowledge, has the aic 
information theoretic been employed to answer the question: How many elements 
yield the ‘best’ Maxwell model for a given data set? 

However, this very question has been asked, and answered, from the viewpoint 
of statistics, where the meter stick has been the minimum of some objective func- 
tion. The outcome of this process is the 7 to 10 mm elements that are typically 
employed, with the actual number of Maxwell chains needed in any given instance 
being dependent upon the actual data being fit. 

We now answer this same question using aic as the meter stick. Aic is a mar- 
riage between statistics and information theory — see App. A for an overview — that 
enables multi-model inference. Presented in Table 3 are the maximum likelihood 
estimates for the parameters in three gmm models with increasing numbers of mm 
elements. Table 4 presents their associated aic statistics. If one were to use the ob- 
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mm elements 

Moo (kPa) 

Mo ( kPa ) 

Cl 

H (s) 

c 2 

t 2 (s) 

C3 

t 3 (s) 

N = 1 

1. 12 

3-45 

1 

1. 16 

— 

— 

— 

— 

N = 2 

0.992 

3-76 

1 -c 2 

0.50 

O.229 

10. 0 

VTjy 

— 

N= 3 

0.861 

3-76 

1 -C2-C3 

0.45 

O.181 

5-05 

0-119 

47.6 


Table 3 Optimized parameters /too.- /-to, T i, C2, T2, C3 and 7:3 for modeling stress relax- 
ation in the human calcaneal fat pad using the gmm kernel function. 


mm elements 

0 

O 

Maic 

A; 

N = 1 

O.II748 

O.IO33 

-8.9 

44-3 

N = 2 

0.00233 

O.OI46 

- 33-7 

! 9-5 

N = 3 

0.0008l 

0.0086 

+9.6 

62.8 


Table 4 Aic statistics for parameter estimates listed in Table 3. 


jective function 0 as the meter stick, or equivalently, the coefficient of variation cr, 
then N = 3 mm units is obviously ‘best’, and it is better than any of the models 
presented in Table 2 . Most likely, this could be improved upon still further by using 
even more mm elements. However, the aic measure for multi-model inference /z AIC 
overwhelmingly selects N = 2 as being the ‘best’ gmm model for the calcaneal fat 
pad. The parameters of the N = 2 gmm model best represent the ‘information’ 
present within the data amongst the various gmm models. Interestingly, this is the 
number of Maxwell chains used by Miller- Young ( 2003 ), where she reported values 
of Tj = 0.5 s and T 2 = 24 s. We are in agreement on the former value but differ on 
the latter. Our differing values for T 2 are likely due to the fact that we obtained our 
parameters from maximum log-likelihood estimates; whereas, Miller- Young ob- 
tained hers from nonlinear regression estimates. We also employed different elastic 
models. Furthermore, the ramp time t\ that she imposed in her analysis was not 
documented. 

Comparing the best gmm model against any of the previous four models via 
their aic differences Aj ranks the best gmm as being a ‘poor’ material model for 
inference according to the criteria put forth in App. A. 


6 Summary 

An elastic strain-energy function has been proposed that has great potential for the 
field of tissue mechanics. An application of the aic information theoretic lead to a 
power-law form of this free energy as being the best choice for the purpose of de- 
scribing compression in the human calcaneal fat pad. The elastic material behavior 
associated with this free-energy function was then analytically continued into the 
thermodynamically irreversible domain of viscoelasticity via the k-bkz hypothesis. 
Because the original k-bkz theory was derived for fluids, and our application is 
for solids, we took the 1 D standard fov solid and converted it into an equivalent 
memory-function format. This gave us the mathematical structure of an elastic/ 
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viscoelastic constitutive equation where strain, not strain rate, is the forcing func- 
tion in the integrand, in accordance with the k-bkz hypothesis. With this overall 
mathematical structure in hand, and with the 3D tensorial structure that the k-bkz 
hypothesis provides (as applied to our elastic strain-energy function), a new class of 
viscoelastic materials was derived. A second application of the aic information the- 
oretic selected the rfd (regularized fractional derivative) as being the best choice 
for the relaxation/ memory function kernels present in our material model for the 
purpose of describing stress relaxation in the fat pads of our feet. 

We have found the aic information theoretic to be a technology of great utility 
in biomechanics applications, yet it is apparently an unknown technology to this 
discipline. It is therefore our hope that biomechanicians will find our explanation of 
it to be straightforward and easy to exploit. Aic provides a means whereby we can 
enhance our understanding of the mathematical models that we use to describe 
the various behaviors that tissues exhibit. 
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A Akaike Information Criterion 


“Truth in the biological sciences and medicine is extremely complicated, and we 
cannot hope to find exact truth or full reality from the analysis of a finite amount of 
data. Thus, inference about truth must be based on a good approximating model. 
Likelihood and least squares methods provide a rigorous inference theory if the 
model structure is ‘given.’ However, in practical scientific problems, the model is 
not ‘given.’ Thus, the critical issue is, ‘what is the best model to use.’ This is the 
model selection problem.” Burnham and Anderson (2002, pg. 47). 

We have used theory to provide mathematical (tensorial) structure to a class of 
material models that contains a known (finite) set of candidates. However, theory 
is unable, at least in our case, to discern which candidate model is ‘best’, especially 
since our models are nonlinear. We therefore desire a methodology whose outcome 
will objectively select the best model from this set of candidate models when fit against 
known data prone to noise. We refrain from subjectively assigning the model, which 
is accepted practice in the biomechanics literature of today. Instead, we employ 
the Akaike Information Criterion (aic) — a technology for use in model selection 
via multi-model inference. Other criteria also exist, see Burnham and Anderson 
(2002, pp. 65—70). Aic is based on the principle of parsimony: A compromise be- 
tween bias-squared (simplicity - increases with decreasing numbers of model pa- 
rameters) and variance (complexity - increases with increasing numbers of model 
parameters). Aic uses maximum log-likelihood inference to obtain ‘optimum’ pa- 
rameter estimates for each candidate model. These estimates, in conjunction with 
the objective function, are then inputs into a Kullback-Leibler (kl) information- 
theoretic that is used to discern the ‘best’ model for inference, selected from the set 
of fitted models. The selected ‘best model’ need not be the ‘model that fits best’. 
Consider an optimization problem where: 

- K is the number of candidate models, 

— L is the dimension of unknown parameters p = \p \ , P2, ■ ■ ■ , 7Fl] t , 

— M is the dimension of state variables y = [v 1 , >’2, . . . , }’m\ , and 

- N is the number of observed variables j 1 = [i’j , 'Tm] T > {*«> }’J 

with tj being the associated times of observation. 

Consider the special case where: 

1. errors between observations y‘ & _v ! + 1 are independent V/ €{ 1 ,...,#—!}, 

2. errors in observations y l are normally distributed about the solution y(tj, p), 
with p being the optimum parameters, 

3. errors between and yl are independent for all k ^ i over all the i, and 

4. a constant coefficient of variation exists in the observations yj , which is inde- 
pendent of j over all the i . 

If the above conditions hold, then Baker et al. (2005) have shown that the maximum 
log-likelihood estimate reduces to a weighted least-squares estimate whose weights 
are elements from the inverse of the covariance matrix of errors, which permits a 
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dimensionless objective function to be defined as 


i=l 7 = 1 V 


y 


(Ai) 


implying a least-squares coefficient of variation of O' = 1 ; whereas, the maximum 
likelihood estimate for the coefficient of variation in the data is given by 




1 

MN 


0(P). 


(Aa) 


Akaike’s (cf. Burnham and Anderson 2002, pp. 60—64) measure for multi-model 
inference is then quantified via 

/A.c = MN In (<P(p)) + 2 (L + 1 ) + ' (As) 

wherein the ‘Pip) of Eq. (Ai) has been minimized to get the maximum likelihood 
estimates p for the model parameters, whose dimension L may vary from model 
to model; however, dimensions M and N remain fixed. The last two terms on the 
right-hand side of /x AIC correct for model bias in the sense of kl information theory. 
The ‘best’ model for the purpose of inference is the one with the smallest or most 
negative p, AIC . 

Confidence intervals can be assigned to each parameter pi in p. If we de- 
note pi = [pi, P2, ■ ■ ■ , Pl -\ . Pi, Pl+\* ■ ■ ■ , Pl] T such that pi € [pf m , p™ x ](x f)> 
then confidence intervals are obtained via the formula (Venzon and Moolgavkar 
1988) 

MN\ln( 0 (p t )) - In ( 0 (p)) \ < x ?, (A4) 

wherein is the / 2 -distribution for 1 degree of freedom, which for the 0.95 quan- 
tile is 3.841, for example. 0 (pi) varies only parameter pi from optimum p in a 
search for those values p™ m and that will satisfy the equality in Eq. (A4). 

For a given data set, a ‘best’ model can be obtained by employing the straight- 
forward methodology outlined above. But will this model be the ‘best’ for another 
data set? Maybe not. Rules have been developed that allow one to dismiss those 
models that are not likely to ever be ‘best’, while retaining a subset of ‘good’ mod- 
els. Begin by constructing the aic differences 

K 

Ai = /x AIC( . - min p AICk . (A5) 

k = 1 

One then applies the following rule to infer which models are ‘good’, which ones 
are ‘mediocre’ and which ones are ‘poor’ (Burnham and Anderson 2002, pg. 70): 


Ai 

Level of Empirical Support for Model i 

0-2 

good model 

4-7 

mediocre model 

> 10 

poor model 
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It is not the absolute size of the aic measure /t. AIC that matters, but rather, it is the 
relative value of the aic difference Aj that is important. The above rule is based on 
the weight of evidence in favor of model i being the actual kl ‘best’ model for the 
problem at hand, given that one of the candidate models is actually this model; in 
other words, this rule has a solid footing in information theory. 


B Other Viscoelastic Kernels 

Besides the fov kernel (a relaxation/ memory function pair) presented in §4, four 
additional viscoelastic kernels are presented in this appendix that complete the set 
of candidate viscoelastic kernels used for selecting the best model for describing 
the dynamic behavior of the human heel pad using the aic methodology presented 
in App. A. These are among the more popular relaxation/ memory function pairs 
that appear in the viscoelastic literature. 


B.i GMM Kernel 


The eminently popular Maxwell model (mm) has a generalized relaxation function 
of a decaying exponential 

G(t) = exp(— t/r), (Bi) 

whose memory function is simply 


M{t) = 


exp(— t/x) 

T 


(Ba) 


with the material constant r (> 0) being called the characteristic time. 

The generalized Maxwell model (gmm) is composed of a finite sum of N dis- 
crete mm elements such that 

N N 

G(t) = ^2 c„ exp(— ?/ r„), 7> = 1, 0 < n < r 2 < • • • < x N , (B 3 ) 

n = 1 n = 1 

whose memory function is 

N 

M(t) = exp i-t/tn), (B4) 

L ' L n 
n= 1 

where each term in the sum can be thought of as being associated with a separate 
integral. It is not always possible to obtain an unique set of parameters for this 
model. The sum over all c n equaling 1 enforces Go = 1 , while G 0 o = 0 follows 
if t„ >0 for all n; hence, the model obeys the principle of fading memory under 
these pretenses. 

Without exception (to our knowledge), gmm is the viscoelastic kernel prepro- 
grammed into commercial finite element codes that have viscoelastic material mod- 
els in them. Gmm is the kernel that arises from a system of first-order differential 
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equations describing viscoelasticity when derived from the theory of internal-state 
variables with N internal variables (cf. Simo and Hughes 1998, Chp. 10). 

Any continuous, linear, viscoelastic spectrum can be discretized, and in doing 
so, can be represented with an approximating gmm kernel. Fulchiron et al. (1993) 
and Simhambhada and Leonov (1993) propose using a Pade-Laplace technique to 
achieve this objective. Here optimum parameters to a Pade expansion of chosen 
order are acquired in the Laplace domain, where the problem is well posed. The 
results are then transformed back into the time domain for use. One usually needs 
about 10 Maxwell chains, i.e., 10 mm kernels or N = 10 , in order to obtain an 
approximation of reasonable accuracy for a continuous spectrum whose frequency 
range is known over 7—10 decades. 

In the thesis of Adolfsson (2003, paper 1), the Voigt fov relaxation spectrum 
was discretized to obtain analytic formulae for the Maxwell chain coefficients C n 
given that x n = nr / N, n = 1,2 , . . . , N , with r being the characteristic relaxation 
time from the Voigt fov model. For a typical value for a of 0.67, he found that 
the normalized relaxation function predicted by 10,000 mm elements to be in about 
1% error with that of the Voigt fov relaxation function, the relaxation function 
obtained by using 1000 mm elements was in about 2% error, and when 100 mm 
elements were used it was in about 5% error. 


B.2 KWW Kernel 


A popular relaxation function from the viscoelastic liquids literature is the stretched 
exponential of Kohlrausch (1847) and Williams and Watts (1970) (kww), which for 
a solid is given by 

G(t) = exp (-(t/vf), (B 5 ) 

whose memory function is 


M(t ) = 


P ex p (-0A/) 


(B6) 


where r (> 0) and P (0 < ft < 1) are the material constants. 

This relaxation function is normalized in the sense that Go = 1 and Goo = 0 . 
The memory function is singular at the origin, i.e., Mq = 00 (given 0 < < 1 ), 

with M{t ) monotonically asymptoting towards Moo = 0 with increasing t . How- 
ever, if /} were allowed to be greater than 1, then Mo = M( 00) = 0 and the 
memory function would no longer be monotonic, violating the principle of fading 
memory. Consequently, 0 < P < 1 in order for Eqs. (B5 & B6) to be in accordance 
with this physical principle. 


B.j QLV Kernel 

Quasi-linear viscoelasticity (qlv) was introduced by Fung (1971), with its relaxation 
function not appearing until much later (Fung 1993, pg. 285). When written as a 
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generalized relaxation function, it becomes 

r-us E l (t/r 2 )- Ei(t/n) 

G(n = — HifeTF) — ’ (I 

with parameters Ti (> 0) and x 2 (> T 1) designating material constants, wherein 


?n(x) := £ 


y~ n e~ xy dy 


is the exponential integral. The qlv relaxation function satisfies G 0 = 1 and 
Goo = 0 . The memory function associated with this relaxation function is more 
user friendly, it being simply 

_ ex P (- ? / r 2) - ex P H/*0 


This has become the de facto standard for characterizing soft-tissue viscoelasticity in 
the biomechanics literature. 

One needs to be careful to distinguish between the Mittag-Leffler function 
E m , n (x ) (especially the one-parameter Mittag-Leffler function E n (x)) and the ex- 
ponential integral E n (x). all of which are their accepted notations. 

The qlv relaxation function is not usually written in the above format. Specifi- 
cally, E 0 o does not appear in the qlv literature; rather, a parameter c (> 0 ) appears 
that relates to the rubbery modulus via E 00 = Eq/[\ + c \n(x 2 /r\)\, where c rep- 
resents the height of a box relaxation spectrum that begins at time Ti and ends at 
time x 2 . Because Mq = \/x\ — \/x 2 'is positive, with M(t ) monotonically decreas- 
ing to o as t — > 00, the Qyv kernel is found to be in accordance with the principle 
of fading memory. 

The paper of Puso and Weiss (1998) employed the gmm model, using 7 mm 
kernels (Eq. B3) to represent the qlv kernel (Eq. B7), so that they could approximate 
£lv in a finite element code in an efficient manner. 

More recently, Doehring et al. (2004) applied the qlv and fov kernels to stress 
relaxation and cyclic data obtained from heart-valve tissues, and found their er- 
rors in predictive capability to be similar. Fov had an advantage over qlv in their 
parameter estimation, as only two of £)lv’s three parameters were observed to be 
sensitive to the data. Parameter x 2 was found to be insensitive, at least to relaxation 
data. This is a well-known fault of qlv. However, we did not experience this diffi- 
culty when fitting the relaxation data for the heel pad, as X 2 was found to lie within 
the time interval of the relaxation experiment. 


/j.y RED Kernel 

Single-integral finite-strain (sifs) viscoelasticity (Johnson et al. 1996) employs a 
relaxation function of the type G(t) = 8/(8 + t) that can be analytically continued 
as a power law so that the relaxation function becomes 


G(t) = 


8 

8 H~ t 


OL 


(Bio) 
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whose affiliated memory function is just 


M{t) = 


ce8 a 

(8 + 0“ +1 ’ 


(Bn) 


where a (> 0 ) and S (> 0 ) are the material constants. Williams (1964) used this 
kernel to describe the relaxation behavior of solid rocket propellants and called it 
the modified power law. 

Here Go = 1 and Goo = 0 , as required, and Mq = a / 8 with M(t) rnono- 
tonically decreasing toward M 0 0 = 0 , thereby ensuring that rfd possesses a fading 
memory kernel. 

This kernel is not an Abel kernel, although it is similar in many respects. Specif- 
ically, Eq. (Bio) is not the Voigt fov kernel G(t) = t~ a / T( 1 — a) associated with 
the fractional derivative in Eq. (19). In the Voigt fov kernel Go = OO and G 0 0 = 0 , 
and therefore, it is singular at the upper limit of integration. Rather, Eq. (Bio) is 
a kind of regularized, fractional derivative (rfd) kernel whose relaxation function G is 
normalized so that Go = 1 . This is accomplished by pushing the singularity just 
outside the interval of integration by a small distance 8 , i.e., the singularity is moved 
to t + 8 . The Voigt fov kernel and the rfd kernel are indistinguishable at large t . 
It is only when t < 8 that these two kernels differ significantly. Exponent a can 
therefore be interpreted as a fractional order of evolution; it is the slope of the re- 
laxation curve through the transition region between glassy and rubbery behavior. 
Similarities and differences between the Voigt fov and rfd kernels have been quite 
thoroughly investigated by Bagley (1987). 

The rfd kernel is not the only way in which a fractional derivative can be 
regularized. Two alternative methods have been proposed in the mathematical lit- 
erature. The first one, completely different from the modified power law of rfd, is 
based on a discretization of the fractional derivative — see, e.g., Tuan and Goren- 
flo (1994a, b). The second one, described by, e.g., Rubin (1996, §11) or Gorenflo 
and Rubin (1994), is much closer to, but not identical with the rfd concept. Their 
method to tackle the singularity in the Voigt kernel mentioned above is very simple: 
Instead of using the full (and singular) integration range from 1 0 to t in the defi- 
nition of the Caputo derivative, Eq. (19), they only integrate from to to t — 8 with 
a certain (positive but small) regularization parameter 8 , thus cutting off the part 
of the interval where the singularity appears; it still occurs at time t. Compared 
to our approach, their method has the charm that the correspondence between 
the kernel ( t — s)~ a and the forcing function f(s) remains unchanged; whereas, 
our scheme shifts one factor by an amount of 6, but does not shift the other factor 
simultaneously. This feature seems to be an advantage of the method of Gorenflo 
and Rubin. On the other hand, their cut-off strategy means that in an actual com- 
putation of the fractional derivative, which is supposed to be a functional with full 
but fading memory, the contribution that is associated with the most recent past 
(the time interval from t — 8 to the current time t) is ignored completely; whereas, 
our method retains this information. 
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C Derivation of Constitutive Model 

Convected tensor fields are employed in this appendix to derive our constitutive 
formulas, which are then transformed into Cartesian tensor equations using well- 
established mapping techniques outlined in §C.2 of this appendix. This process of 
deriving constitutive expressions using convected tensor fields, and then mapping 
them into Cartesian space, ensures that the resulting Cartesian formulae are frame 
invariant. The interested reader is referred to the classic paper by Oldroyd (1950) 
and the various textbooks of Lodge (1964, 1974, 1999) for a thorough treatment of 
this subject matter. Apparently, Zaremba (1903) was the first to realize that the use 
of convected fields automatically ensures frame indifference. 


C.i Constitutive Formulation 


Consider a mass element £ whose density is Q = q(% ; t), which is located by a set 
of components (£* , £ 2 , f 3 ) in a convected coordinate system whose numeric values 
do not vary over time. Mass point £ is put into a state of stress jr ,J = 7 t ,J ((• ; t) as 
the result of a change in shape d y,-j := yij (£ ; t + df ) — yij (§ ; t ) imposed over an 
increment in time df. This deformation induces a differential change in the work 
done d W = d W(%', t,dt) on the mass element (including an energetic contribu- 
tion associated with the kinetic energy of the mass point) that is quantified by the 
formula (Lodge 1974, pp. 194— 195) 

dW = i- n ij d yji , (Ci) 

where repeated indices — one a subscript, the other a superscript — are summed 
over in the usual manner. The stress tensor n has contravariant components 7 t ,J , 
the metric tensor y has covariant components yij, while the inverse metric y ~ 1 
has contravariant components y lJ . y~ l exists because y is symmetric and positive 
definite by definition, viz., d.y 2 (f; ; t ) = df ' (§ ) y/j (§ ; t ) df 7 (§ ) where d.v (> 0) is 
the length of infinitesimal vector d§ at particle § and time t. 

The metric tensor t) of convected tensor analysis is a function of time. 
This characteristic is not shared by the metric tensor g (jc ) of general tensor analy- 
sis, which is independent of time (c£, e.g, with Sokolnikoff 1964). 

Unlike Cartesian fields, which are evaluated in fixed, rectangular-Cartesian, 
coordinate systems, convected body fields are evaluated in curvilinear coordinate 
systems that move (i.e., convect) with the body as though the coordinate axes were 
drawn onto the body. 

Tensors JC and y are the fundamental fields of convected tensor analysis from 
which constitutive equations are constructed. 

From thermostatics, Eq. (Ci) leads to a potential-based constitutive equation 
for the elastic state of stress characterized by (Lodge 1964, pp. 154-161) 


f n ij 


f dW 3W \ 

V fyij d'/ji ) ’ 


TC ‘ 


TC 


Ji 


(C2) 
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where Q o = p(§; to). From the conservation of mass, one is lead to the expression 
TT = ( det (I'd V)) V2 wherein the tensor contraction y Q 1 y has mixed compo- 
nents y l/c ( £ ; to)ykj(t ; ; t). Symmetry in stress, a consequence of the thermostatic 
potential W, is also in accordance with the conservation of angular momentum. 

Following the suggestion of Flory (1961), a hydrostatic/ deviatoric split in the de- 
formation held is introduced in such a way that the strain energy can be decoupled 
as W = Wj, + Wj, where W b and Wj are the bulk and deviatoric strain energies, 
respectively. This allows us to replace Eq. (C2) with a decomposition in the state of 
stress given by 

71 ij + p y ij = n ij , n iJ = n ji , n kt ytk = o, (c 3 ) 


wherein the hydrostatic pressure p := — 
equation 


P = “Co 


j 7 t kt yik is governed by the constitutive 


dWb(J) 

dj 


(C 4 ) 


with J := (det()/Q 1 y)) ~ quantifying the volume of mass point £ , in a relative 
sense. The extra stress II has contravariant components fl ,j that are deviatoric 
(i.e., rjk^yik = 0), and is governed by the constitutive equation 


f n ij 


Qo J 2 ^ Dev 


'dW d (yo,y) + dW d (y 0 .y)' 


dYij 


d Yji 


(C 5 ) 


with Dev [•]' / := (•) ,y — y [{■) k( 'Yik\ Y‘ J dehning the deviatoric components of a 
contravariant held. (We will not need the deviatoric operators for covariant and 
mixed tensor helds.) Tensor y has components py := J 2 ^ yij (§ ; t ), and is there- 
fore isochoric because det(yg*p) = 1 . The argument y 0 in Wj is a constant 
tensor held establishing the initial shape of the mass element. It is like a built-in 
boundary condition. 

Soft connective tissues are inherently viscoelastic. By adopting the design phi- 
losophy advocated in the k-bkz theory of viscoelasticity (Kaye 1962; Bernstein 
et al. 1963), one can analytically continue the above elastic state into a viscoelastic 
state by employing an appropriate expression for the strain-energy gradient as the 
forcing function in the viscoelastic structure of Boltzmann (1874). 

In soft-tissue mechanics, it is reasonable to assume that only the deviatoric re- 
sponse is viscoelastic. There are applications where viscoelastic compressibility can 
be very important (cf. Leonov 1996); however, the in vivo rate-controlling relaxation 
mechanisms of soft tissues are not known to be affiliated with volume change. 


C.1.1 Bulk behavior. Soft tissues are nearly incompressible materials in that their 
bulk moduli are orders of magnitude greater than their shear moduli, and as such, 
it is reasonable to consider a convex pressure/volume model whose spherical en- 
ergy is given by (Sirno and Hughes 1998, pg. 361) 

eo W b {J) = K\(\(j 2 -l)-\nj), 


(C6) 
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which, from Eq. (C4), yields a symmetric expression for hydrostatic pressure of the 
form 

P = ~ K W~ /_1 )’ ( c 7) 

with k being the bulk modulus. Dilatation j(J — / _1 ) is a second-order accurate 
approximation to the dilatation of Hencky (1928), viz., In J (cf. Freed 2004). 

Again, the viscoelastic attributes of bulk behavior are not considered to be rate 
controlling in soft-tissue mechanics. Bulk behavior is taken to be elastic, and in 
many cases can be assumed to be incompressible. 


C.1.2 Deviatoric behavior. Using Eq. (27) as our template (i.e., the memory function 
formulation of ID fov), and using the k-bkz hypothesis to establish mathematical 
structure, we postulate the existence of the following general constitutive equation 
for the deviatoric response of a class of viscoelastic solids 


f n ij = (/too + (/to - /too) G(t)) 
x Qo J~ 2 ^ Dev 


dW d (y 0 ,y) + dW d (y 0 ,y) 


f 

J t(\ 


9 Yij 

M(t—s ) q s j ~ 2 ' 3 Dev 


9 Yji 


+ (MO - /too) 


dW d (y s ,y) + dW d (y s ,y) 


9 Yij 


9 Yji 


d.v, (C8) 


where /to and /too are the glassy and rubbery shear-moduli, respectively, with 
y := J~ 2 ^y and y := J~ 2 ^y, where J := (det(y“ 1 y)) ' 2 with y s having com- 
ponents yij(^',s). Metrics y and y are isochoric in the sense that det(j/Q*y) = 
detfy” 1 y) = 1 ; furthermore, y = y whenever s = to- The extra stress FI is com- 
pletely described here by two material functions: W d and G, because M(t—s ) = 
dG(t—s)/ds by definition. 

The integral for extra stress in the above expression is in accordance with the 
k-bkz hypothesis: The affect of configuration s < t on stress IJ at time t is equiv- 
alent to the affect of elastic energy stored by the material with the configuration at 
time s serving as its reference state, weighted by a memory function and summed 
over the history of all past configurations. The integral in the theory of Kaye (1962) 
and Bernstein et al. (1963) goes from —00 to t ; whereas, our integral goes from 1 0 to 
t , with time to being a stress-free equilibrium state. Their theory is for viscoelastic 
fluids; our theory is for viscoelastic solids. This is why our formulation has an elas- 
tic non-integral contribution, and why ours has an initial reference configuration 
associated with time to . 

The constitutive assumption that we employ to describe the isotropic deviatoric 
response of biological tissue is Eq. (13). When rewritten for an arbitrary reference 
state s € [to. t\, and in terms of converted fields, it becomes 


2£>o W d (y 0 .y) = p,j(f(pGl)'*rf(pi;3) + f{p 2 \U)-f(p 2 \ 3)), 

(C9) 
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whose invariants are defined by 

^ = tr(^o 'y)- # = tr(y“V 0 ), 

A. ^ ^ (Cio) 

/ = tr(y;^), ff = a{f- 1 y s ). 

We observe in Eqs. (Cg & Cio) that the assigned strain-energy function weighs the 
two states S and t equally, within a constant of proportionality introduced via p \ 
and p2- Function f is subject to the stability criteria listed in Eq. (14). After some 
algebra, the above definitions transform Eq. (C8) into the following constitutive 
model 


r n = 2(/too + (Mo - Moo) G(t)) 

x j~ 2/i Dev [\{f'(p\ \ 1 ) y^ ] - f'{P 2 \n)y~ l y 0 y~')] + 2(mo-Moo) 
x J M(t-s) J- 2/ i Dev [\(f( P ul) y; ] - f(p 2 \ ff) p~' y s p~') 


d.v. 

(C11) 


We leave the functional forms for f and G as unspecified at this time, recalling 
that M(t — s) = 3 G(t — s)/ds. 


C.2 Field Transfer 


Field transfer is a process by which convected (body) tensor fields are mapped into 
general (space) tensor fields or, as in our case, into Cartesian tensor fields. This 
has been thoroughly documented by Lodge (1964, 1972, 1974), and summarized by 
Freed (1995). 

Because the mapping of a convected body field into a Cartesian space field is 
many-to-one, we use the notation =r- to denote it. The underlying mathematics of 
this operation are substantial and are not detailed here. They can be found in the 
references cited above for the reader who is interested. The field transfer operator 
is time dependent. Whenever a mapping is to be done at current time t, denoted as 

=>, the resulting Cartesian fields will be defined in the Eulerian frame. Likewise, 

to 

whenever field transfer is to occur at reference time to, denoted as =>, the ensuing 
Cartesian fields will be defined in the Lagrangian frame. 

Whenever the time dependence of the convected metric tensor y coincides 
with the time of field transfer, the resulting Cartesian fields are the identity tensor I; 
specifically, 


t 





y 0 1 



(C12) 


which illustrates the many-to-one property of this mapping. The covariant, con- 
travariant or mixed nature of a convected tensor field is lost by the field-transfer 
operator whenever a convected body field is mapped into a Cartesian spatial field. 
This property is not lost, however, whenever convected body fields are mapped into 
general spatial fields, which are not used in this paper. 
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In contrast with the mappings of Eq. (C12), deformation fields will result when- 
ever the time dependence of the converted metric differs from the time of held 
transfer; in particular, 

Fo=> b - Yo * = ^' B - Y => C, y 1 => C \ (C 13 ) 


where B and C are the left and right (or Finger 1894 and Green 1841) deformation 
tensors, respectively, defined in Eq. (3). 

Similarly, whenever metric fields are referenced to the dummy state of integra- 
tion s, to < s < t , for example, as they can appear in viscoelastic models, then 
it follows straightaway that these body metric tensors map into Cartesian space 
according to 


Y S ^ B- 1 , y s ^U C, y~ l ^ C\ (C14) 

where B = FF T and C = F T F, with F and F being dehned in Eq. (2). 

The convected stress tensor n maps into Cartesian space as it => T so that 


^-1 


P := 


and n 


1 0 


S := f F 


-1 


TF 


(C15) 


where T is the Cauchy stress, P is the Kirchhoff stress and S is the second Piola- 
Kirchhoff stress, all of which are symmetric. 

Applying the field-transfer results given in Eqs. (C12-C15) to the constitutive 
formulas listed in Eqs. (C3, C7 & Cii) produces the constitutive formulas written 
down in Eqs. (30—34). 


D Small Strain/Rotation Relationships 

A variant of one’s theory is often sought that can be expressed in terms of the 
displacement vector 1/ := x (A ) — X and its gradient G := F — I , G,y = <)Uj / 3 Xj , 
through the classic tensor fields defining infinitesimal strain e := \{G + G t ) and 
infinitesimal rotation co := ^-(G — G T ), wherein the strain tensor e is symmetric 
(i.e., e = f T ), while the rotation tensor co is skew (viz., co = — « T ). 

From these definitions, one can express the finite stretch and rotation tensors 
U and R (where F = RU from the polar decomposition theorem) in terms of the 
infinitesimal strain and rotation tensors e and co via the formulae 

U = exp (e) = I + t + \e 2 + ■ ■ ■ , 

U -1 = exp (— f ) = I-( + jf 2 , 

R = exp (o>) = I + co + jco 2 -I , 

R t = exp (—co) = I — co + jco 2 , 


(Di) 
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so that the various forms of the deformation gradient become 

F = RU = I + t + « + &) t + i(e 2 + w 2 ) H , 

F t = UR t = I + f - a - e a + \{e 2 + co 2 ) + ■ ■ ■ , 

F _1 = U _ 1 R t = l - e - m + e m + j(e 2 + co 2 ) , 

F- t = RU -1 = I - t + at - co f + \(e 2 + co 2 ) , 

which in turn produce the following expressions for the deformation tensors 


C =F T F = I + 2e + 2e 2 + -- - , 

C -1 = F"'F“ t = I — 2e + 2t 2 , 

B = FF T = I + 2e + 2(o> t — e co) + 2e 2 + ■ ■ ■ , 

B " 1 = F _T F _1 = I — 2e — 2(co e - e to) + 2e 2 . 


Using these series expansions allows one to construct first- and second-order ap- 
proximations for most finite-strain constitutive formulae. 
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